a <- read.csv("test.csv",row.names=1,header=T)
sample <- as.factor(c(rep("Control", 3), rep("one", 3), rep("three", 3), rep("five", 3)))
res.Control <- selectHKgenes(a[sample=="Control",], method = "Vandesompele", geneSymbol = names(a), minNrHK = 2, trace = T, na.rm = T)
res.one <- selectHKgenes(a[sample == "one",], method = "Vandesompele", geneSymbol = names(a), minNrHK = 2, trace = TRUE, na.rm = TRUE)
res.three <- selectHKgenes(a[sample == "three",], method = "Vandesompele", geneSymbol = names(a), minNrHK = 2, trace = T, na.rm = TRUE)
res.five <- selectHKgenes(a[sample == "five",], method = "Vandesompele", geneSymbol = names(a), minNrHK = 2, trace = FALSE, na.rm = TRUE)
ranks <- data.frame(c(1, 2:5), res.Control$ranking, res.one$ranking, res.three$ranking, res.five$ranking)
names(ranks) <- c("rank", "Control", "1 mM", "3 mM", "5 mM")
print(ranks)

#install RColorBrewer from CRAN > require(RColorBrewer) 
#stability of gene expression upon exclusion of less stable genes
mypalette <- brewer.pal(4, "Set1")
matplot(cbind(res.Control$meanM, res.one$meanM, res.three$meanM, res.five$meanM), type = "b", ylab = "Average expression stability M", xlab = "Number of remaining control genes", axes = FALSE, pch = 19, col = mypalette, ylim = c(0.005, 0.045), lty = 1, lwd = 2, main = "Gene stability 6h")
axis(1, at = 1:4, labels = as.character(5:2))
axis(2, at = seq(0.005, 0.045, by = 0.01), labels = as.character(seq(0.005, 0.045, by = 0.01)))
box()
abline(h = seq(0.005, 0.035, by = 0.01), lty = 2, lwd = 1, col = "grey")
legend("topright", legend = c("No treatment","1 mM NaF","3 mM NaF","5 mM NaF"), fill = mypalette)

#variation of expression upon inclusion of additional genes beyond 3 in the calculation
mypalette <- brewer.pal(3, "YlGnBu")
barplot(cbind(res.Control$variation, res.one$variation, res.three$variation, res.five$variation), ylab="Pairwise variation",ylim=c(0,0.0125),beside = TRUE, col = mypalette, space = c(0, 2), names.arg = c("No treatment", "1 mM NaF", "3 mM NaF", "5 mM NaF"),main="Variation of gene expression 6h")
legend("topright", c("V4/5", "V3/4", "V2/3"), xpd = TRUE, horiz = F, inset = c(0, 0), fill = mypalette, ncol=3)
abline(h = seq(0, 0.01, by = 0.003), lty = 2, col = "grey")
